/*!  \file  RandomGeneraor.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      05/12/2013 06:09:56 PM
 *   \copyright GNU Public License.
 */
#include <assert.h>
#include <math.h>
#include <string.h>
#include <stdint.h>
#include <stdlib.h>
#include <iostream>
using namespace std;
#include "RandomGenerator.h"
RandomGenerator::~RandomGenerator(){ assert(stateData == 0);}
void RandomGenerator::setSeed(const unsigned int N) {
	// stateData pointer should always be alligned with int
	assert(reinterpret_cast<uintptr_t>(stateData) % sizeof(N) == 0);
	size_t s = stateSize();
	int limit = s / sizeof(N); 
	limit += s % sizeof(N) != 0;
	for (int i = 0; i < limit; i++)
		reinterpret_cast<unsigned int *>(stateData)[i] = N;
}
void RandomGenerator::setSeed(void * stateDataPtr) {
#ifndef NDEBUG
	void * res = 
#endif
	memcpy(stateData, stateDataPtr, stateSize());
	assert(res == stateData);
}
void RandomGenerator::fillUniform(const unsigned int N, double * __restrict__ d) {
	for (unsigned int i = 0; i < N; i++)
		d[i] = getUniform();
}

// drand48 family of generators
class Rand48 : public RandomGenerator {
public:
	Rand48(const unsigned int someSeed) {
		size_t sizeInMem = stateSizeInMem();
		assert (sizeInMem % sizeof(int) == 0); 
		unsigned int * p = new unsigned int[sizeInMem /sizeof(unsigned int)];
		stateData = reinterpret_cast<void *>( p);
		if (stateData) setSeed(someSeed);
	}
	virtual ~Rand48() { 
		unsigned int * p = reinterpret_cast<unsigned int *>(stateData);
		stateData = 0;
		if (p) delete [] p; p = 0;
	}
	virtual size_t stateSize() const {return 3 * sizeof(unsigned short int);};
	double getUniform(); // returns random number from in [0.0,1.0)
};
double Rand48::getUniform(){
	unsigned short * state =  reinterpret_cast<unsigned short *>(stateData); 
	unsigned int res = nrand48(state);
	const double rm = 1.0 / ( 256 * 256 * 256) * (1.0/ 128);
	return rm * res;
};
// implements Zaman & Marsaglia ranmar generator
class Ranmar : public RandomGenerator {
public:
	Ranmar(const unsigned int someSeed = 1802) {
		stateData = 0; 
		setSeed(someSeed);
//		high_res = false;
	}
	virtual ~Ranmar() {};
	double getUniform(); // returns random number from in [0.0,1.0)
	virtual size_t stateSize() const { return 102 * sizeof(int);};
	void fillUniform(const unsigned int N, double * __restrict__ d);
	void setSeed(const unsigned int N) { setState(N, 9373);};
private:
	void setState(int ixx, int jxx);
   int        iseed1, iseed2;
   int        ix, jx, c;
   int        u[97];
   static     int cd, cm;
   static     double  twom24;
//   bool       high_res;
};
double Ranmar::twom24 = 1./16777216.;
int Ranmar::cd = 7654321;
int Ranmar::cm = 16777213;
void Ranmar::setState(int ixx, int jxx) {
	if( ixx <= 0 || ixx >= 31328) ixx = 1802;
	if( jxx <= 0 || jxx >= 30081) jxx = 9373;
	iseed1 = ixx; iseed2 = jxx;
	int i = (ixx/177)%177 + 2;
	int j = ixx%177 + 2;
	int k = (jxx/169)%178 + 1;
	int l = jxx%169;
	for(int ii=0; ii<97; ii++) {
   	int s = 0; int t = 8388608;
		for(int jj=0; jj<24; jj++) {
      	int m = (((i*j)%179)*k)%179;
       	i = j; j = k; k = m;
        	l = (53*l+1)%169; 
        	if( (l*m)%64 >= 32 ) s += t;
        	t /= 2;
    	}
		u[ii] = s;
	}
	c = 362436; ix = 96; jx = 32;
}
double Ranmar::getUniform() {
	register int r = u[ix] - u[jx--];
	if( r < 0 ) r += 16777216;
	u[ix--] = r;
	if(ix < 0) ix = 96;
	if(jx < 0) jx = 96;
	c -= cd; if( c < 0 ) c+=cm;
	r -= c; if( r < 0 ) r+=16777216;
	double res = twom24*r;
/* if( high_res ) {
   	register int r = u[ix] - u[jx--];
      if( r < 0 ) r += 16777216;
      u[ix--] = r;
      if(ix < 0) ix = 96;
      if(jx < 0) jx = 96;
      c -= cd; if( c < 0 ) c+=cm;
      r -= c; if( r < 0 ) r+=16777216;
      res += twom24*twom24*r;
	} */
	return res;
}

void Ranmar::fillUniform(const unsigned int n, double * __restrict__ d) {
    for(unsigned int i = 0; i < n; i++) {
        register int r = u[ix] - u[jx--];
        if( r < 0 ) r += 16777216;
        u[ix--] = r;
        if(ix < 0) ix = 96;
        if(jx < 0) jx = 96;
        c -= cd; if( c < 0 ) c+=cm;
        r -= c; if( r < 0 ) r+=16777216;
        d[i] = twom24*r;
    }
/*     if( high_res ) {
        for(unsigned int i=0; i<n; i++) {
            register int r = u[ix] - u[jx--];
            if( r < 0 ) r += 16777216;
            u[ix--] = r;
            if(ix < 0) ix = 96;
            if(jx < 0) jx = 96;
            c -= cd; if( c < 0 ) c+=cm;
            r -= c; if( r < 0 ) r+=16777216;
            d[i] += twom24*twom24*r;
        }
    }
*/
}

RandomGenerator* RandomGenerator::defaultGenerator(const unsigned int someSeed ) {
	return dynamic_cast<RandomGenerator*>(new Ranmar(someSeed));
//	return dynamic_cast<RandomGenerator*>(new Rand48(someSeed));
}
#ifndef NDEBUG
#include "UnitTester.h"
#include <algorithm>
UNITTEST(RandomGeneratorTest);
RandomGeneratorTest::RandomGeneratorTest() {
	setName 				("Test Random Generator");
	setDescription  	("Tests if random generator behaives correctly");
	setTestAuthor		("Yordan Madzhunkov");
	setDeveloperAuthor("Yordan Madzhunkov");
	setCompileTime(__DATE__, __TIME__); // each class that inherits this one should also have this line.
};
#define SQR(x) (x*x)
	double pks(double z) {
		assert(z >= 0.);
		if (z == 0.) return 0.;
		if (z < 1.18) {
			double y = exp(-1.23370055013616983/SQR(z));
			return 2.25675833419102515*sqrt(-log(y))
				*(y + pow(y,9) + pow(y,25) + pow(y,49));
		} else {
			double x = exp(-2.*SQR(z));
			return 1. - 2.*(x - pow(x,4) + pow(x,9));
		}
	}

	double qks(double z) {
		assert(z >= 0.);
		if (z == 0.) return 1.;
		if (z < 1.18) return 1.-pks(z);
		double x = exp(-2.*SQR(z));
		return 2.*(x - pow(x,4) + pow(x,9));
	}

// Small values of prob show that the cumulative distribution function of data
// is significantly different from func.
double KolmogorovSmirnovUniform(const unsigned int N, double * data) {
	int j = 0, n = N;
	double dt,en,ff,fn,fo=0.0;
	sort(data, data + N);
	en=n;
	double d = 0.0;
	for (j=0;j<n;j++) {
		fn=(j+1)/en;
		ff= data[j];
		dt=fmax(fabs(fo-ff),fabs(fn-ff));
		if (dt > d) d=dt;
		fo=fn;
	}
	en=sqrt(en);
	return qks((en+0.12+0.11/en)*d); 
}
bool RandomGeneratorTest::test() const {
	RandomGenerator * rng = RandomGenerator::defaultGenerator(1312);
	const unsigned int N = 12624;
	double * data = new double[N]; if (!data) return false;
	rng->fillUniform(N, data);
	double p = KolmogorovSmirnovUniform(N, data); // Kolmogorov-Smirnov Test
/* 
	cout << endl << p << " for N = " << N << endl;
	unsigned int nBins = 11;
	unsigned int count[nBins]; 
	memset(count, 0, sizeof(count));
	unsigned int bin = 0;
	for (unsigned int i = 0; i < N; i++)
		if (data[i] < (bin+1) * 1.0/ nBins) {
			count[bin] += 1;			
		} else {
			while (data[i] >= (bin+1) * 1.0/ nBins) bin++;
			count[bin] += 1;			
		}	
	assert(bin < nBins);
	for (unsigned int i = 0; i < nBins; i ++)
		cout << count[i] << " ";	
	cout << endl;
*/
	if (data) delete [] data; data = 0;
	delete rng;
	return p > 0.95;
}
RandomGeneratorTest randomGeneratorTest;
#endif
